####Ethnic Bias in Criminal Sentencing - R Analysis File####
###Yue Hou, Assistant Professor, University of Pennsylvania Department of Political Science, yuehou@sas.upenn.edu###
###Rory Truex, Assistant Professor, Princeton University Department of Politics and Woodrow Wilson School of Public and International Affairs, rtruex@princeton.edu###

###LOAD PACKAGES AND FUNCTIONS#
setwd('~/ReplicationFiles')
rm(list=ls(all=TRUE))

library("rvest")
library("cobalt")
library("dplyr")
library("lmtest")
library("sandwich")
library("MASS")
library("qdap")
library("rvest")
library("tidyverse")
library("cobalt")
library("ggplot2")
library("splines")
require("broom") 
library("ids")
library("plyr")
library("grid")
library("readxl")
library("geojsonio")
library("ggplot2")
library("tidyverse")
library("viridis")
library("lme4")
library("RCurl")
library("readr")
library("lfe")
library("psych")

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
library(gridExtra)

require(grid)
vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
arrange_ggplot2 <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
  dots <- list(...)
  n <- length(dots)
  if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
  if(is.null(nrow)) { nrow = ceiling(n/ncol)}
  if(is.null(ncol)) { ncol = ceiling(n/nrow)}
  
  grid.newpage()
  pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
  ii.p <- 1
  for(ii.row in seq(1, nrow)){
    ii.table.row <- ii.row  
    if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
    for(ii.col in seq(1, ncol)){
      ii.table <- ii.p
      if(ii.p > n) break
      print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
      ii.p <- ii.p + 1
    }
  }
}

###PREPARING DATASETS###
load('data.RData')
data<-data2
data$ind<-1
data.year<-aggregate(list(data$ind), by = list(data$ruling.year), sum)
colnames(data.year)<-c("year","cases")


##############################################
#Fig A2: Drug cases in China Judgements online
##############################################
pdf('fig-cases-year.pdf', width=6.23, height=3.5)  
data.year<-subset(data.year, data.year$year!="2018")
ggplot(data.year, aes(x=year, y=cases,color="dodgerblue3"))   + ylab("Total Cases")  + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE) + theme_bw() + theme(legend.title=element_blank()) + scale_colour_manual(values=c("dodgerblue3")) + geom_point() + geom_line()  + theme(legend.position = "none") + scale_x_continuous(name="Ruling Year", breaks=c(2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017), labels=c(2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017), limits=c(2006,2017)) 
dev.off() 

data.filter<-subset(data, data$ruling.year>2013 &  data$def.multiple!=1 & data$multiple_def!=1 & crime.act.aiding!=1 & crime.act.growing!=1 & data$crime.drug!="marijuana" & data$crime.drug!="cocaine" & data$crime.drug!="opium" & data$crime.drug!="other") 

data.filter$def.edu.low<-0
data.filter$def.edu.low[grepl("文盲", data.filter$def.edu)==TRUE]<-1
data.filter$def.edu.low[grepl("小学", data.filter$def.edu)==TRUE]<-1
data.filter$def.edu.low[data.filter$def.edu==""]<-NA
ftable(data.filter$def.edu.low)
mean(data.filter$def.edu.low, na.rm=TRUE)

mean(data.filter$def.mental)
ftable(data.filter$def.mental)

data.meth<-subset(data.filter, data.filter$drug.meth==1)
data.heroin<-subset(data.filter, data.filter$drug.heroin==1)

data.anhui<-subset(data.filter, data.filter$province=="Anhui")
data.beijing<-subset(data.filter, data.filter$province=="Beijing")
data.chongqing<-subset(data.filter, data.filter$province=="Chongqing")
data.fujian<-subset(data.filter, data.filter$province=="Fujian")
data.gansu<-subset(data.filter, data.filter$province=="Gansu")
data.guangdong<-subset(data.filter, data.filter$province=="Guangdong")
data.guangxi<-subset(data.filter, data.filter$province=="Guangxi")
data.guizhou<-subset(data.filter, data.filter$province=="Guizhou")
data.hainan<-subset(data.filter, data.filter$province=="Hainan")
data.hebei<-subset(data.filter, data.filter$province=="Hebei")
data.heilongjiang<-subset(data.filter, data.filter$province=="Heilongjiang")
data.henan<-subset(data.filter, data.filter$province=="Henan")
data.hubei<-subset(data.filter, data.filter$province=="Hubei")
data.hunan<-subset(data.filter, data.filter$province=="Hunan")
data.innermongolia<-subset(data.filter, data.filter$province=="Inner Mongolia")
data.jiangsu<-subset(data.filter, data.filter$province=="Jiangsu")
data.jiangxi<-subset(data.filter, data.filter$province=="Jiangxi")
data.jilin<-subset(data.filter, data.filter$province=="Jilin")
data.liaoning<-subset(data.filter, data.filter$province=="Liaoning")
data.ningxia<-subset(data.filter, data.filter$province=="Ningxia")
data.qinghai<-subset(data.filter, data.filter$province=="Qinghai")
data.shaanxi<-subset(data.filter, data.filter$province=="Shaanxi")
data.shandong<-subset(data.filter, data.filter$province=="Shandong")
data.shanghai<-subset(data.filter, data.filter$province=="Shanghai")
data.shanxi<-subset(data.filter, data.filter$province=="Shanxi")
data.sichuan<-subset(data.filter, data.filter$province=="Sichuan")
data.tianjin<-subset(data.filter, data.filter$province=="Tianjin")
data.tibet<-subset(data.filter, data.filter$province=="Tibet")
data.xinjiang<-subset(data.filter, data.filter$province=="Xinjiang")
data.yunnan<-subset(data.filter, data.filter$province=="Yunnan")
data.zhejiang<-subset(data.filter, data.filter$province=="Zhejiang")

ftable(data.filter$crime.drug, data.filter$crime.act)

#################################
#TABLE A1. Descriptive statistics 
#################################
describe(data.filter[ , c("pun.severity", "pun.lifedeath","def.minority", "def.age","def.female","def.mental","def.edu.low","def.recid","drug.pooled.quantity","international", "minors", "def.pleadnotguilty","def.goodattitude")], fast=TRUE)

##VISUAL ANALYSIS##

######################################################################
#Figure 1: Drug Quantity, Ethnicity, and Sentencing Outcomes in Yunnan
######################################################################

data.yunnan.meth<-subset(data.meth, data.meth$province=="Yunnan")
data.yunnan.heroin<-subset(data.heroin, data.heroin$province=="Yunnan")

data.yunnan.meth.noposs<-subset(data.meth, data.meth$province=="Yunnan")
data.yunnan.heroin.noposs<-subset(data.heroin, data.heroin$province=="Yunnan")

#pdf('fig-quantity-noposs-scatterplot-twopanel-yunnan.pdf', width=6.5, height=8)  
#p1<-ggplot(data.yunnan.meth.noposs, aes(x=drug.meth.quantity, y=pun.severity, color=bin.minority)) + xlab("Drug Quantity - Methamphetamine (grams)") + ylab("Punishment Severity (months)")  + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE) + theme_bw() + theme(legend.title=element_blank()) + scale_colour_manual(values=c("dodgerblue3","dodgerblue3","dodgerblue3","brown3","brown3","brown3")) + annotate("rect", xmin = 0, xmax = 10, ymin = 0, ymax = 84, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 10, xmax = 50, ymin = 84, ymax = 180, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 50, xmax = 200, ymin = 180, ymax = 400, alpha = .15, fill="grey50",color="grey50") + geom_point(size=.15) + stat_smooth(method="loess", se=FALSE) + scale_x_continuous(limits = c(0, 200))+ scale_y_continuous(limits = c(-1, 400)) + geom_vline(xintercept = 10, linetype="dashed", size=.8) + geom_vline(xintercept = 50, linetype="dashed", size=.8, color="black") + theme(legend.position = "none") +  annotate("rect", xmin = 135, xmax = 197, ymin = 1, ymax = 71.5, fill="white", color="grey50") +  annotate("rect", xmin = 139, xmax = 145, ymin = 9, ymax = 19, alpha = .15, fill="grey50", color="grey50") +  annotate("text", label="Recommended sentence", x=170, y=15, size=3.25) + annotate("text", label="Han defendants", x=162, y=36, size=3.25) +  annotate("text", label="Minority defendants", x=165.5, y=56, size=3.25)+ annotate("text", label="Death sentence", x=12, y=375, size=2.75) + annotate("text", label="Life imprisonment", x=14, y=280, size=2.75) + annotate("text", label="15 Year fixed term imprisonment", x=25, y=195, size=2.75) + geom_point(aes(x=142.5, y=56), colour="brown3") + geom_point(aes(x=142.5, y=36), colour="dodgerblue3")
#p2<-ggplot(data.yunnan.heroin.noposs, aes(x=drug.heroin.quantity, y=pun.severity, color=bin.minority)) + xlab("Drug Quantity - Heroin (grams)") + ylab("Punishment Severity (months)")  + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE) + theme_bw() + theme(legend.title=element_blank()) + scale_colour_manual(values=c("dodgerblue3","dodgerblue3","dodgerblue3","brown3","brown3","brown3")) + annotate("rect", xmin = 0, xmax = 10, ymin = 0, ymax = 84, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 10, xmax = 50, ymin = 84, ymax = 180, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 50, xmax = 200, ymin = 180, ymax = 400, alpha = .15, fill="grey50",color="grey50") + geom_point(size=.15) + stat_smooth(method="loess", se=FALSE) + scale_x_continuous(limits = c(0, 200))+ scale_y_continuous(limits = c(-1, 400)) + geom_vline(xintercept = 10, linetype="dashed", size=.8) + geom_vline(xintercept = 50, linetype="dashed", size=.8, color="black") + theme(legend.position = "none") +  annotate("rect", xmin = 135, xmax = 197, ymin = 1, ymax = 71.5, fill="white", color="grey50") +  annotate("rect", xmin = 139, xmax = 145, ymin = 9, ymax = 19, alpha = .15, fill="grey50", color="grey50") +  annotate("text", label="Recommended sentence", x=170, y=15, size=3.25) + annotate("text", label="Han defendants", x=162, y=36, size=3.25) +  annotate("text", label="Minority defendants", x=165.5, y=56, size=3.25)+ annotate("text", label="Death sentence", x=12, y=375, size=2.75) + annotate("text", label="Life imprisonment", x=14, y=280, size=2.75) + annotate("text", label="15 Year fixed term imprisonment", x=25, y=195, size=2.75) + geom_point(aes(x=142.5, y=56), colour="brown3") + geom_point(aes(x=142.5, y=36), colour="dodgerblue3")
#arrange_ggplot2(p1,p2, nrow=2, ncol=1)
#dev.off() 

pdf('fig-quantity-noposs-scatterplot-twopanel-yunnan-loess.pdf', width=6.9, height=8)  
p1<-ggplot(data.yunnan.meth.noposs, aes(x=drug.meth.quantity, y=pun.severity, color=def.minority.fac)) + xlab("Drug Quantity - Methamphetamine (grams)") + ylab("Punishment Severity (months)")  + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE) + theme_bw() + theme(legend.title=element_blank()) + scale_colour_manual(values=c("dodgerblue3","brown3")) + annotate("rect", xmin = 0, xmax = 10, ymin = 0, ymax = 84, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 10, xmax = 50, ymin = 84, ymax = 180, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 50, xmax = 200, ymin = 180, ymax = 400, alpha = .15, fill="grey50",color="grey50") + geom_point(size=.15) + stat_smooth(method="loess", se=FALSE) + scale_x_continuous(limits = c(0, 200))+ scale_y_continuous(limits = c(-1, 400)) + geom_vline(xintercept = 10, linetype="dashed", size=.8) + geom_vline(xintercept = 50, linetype="dashed", size=.8, color="black") + theme(legend.position = "none") +  annotate("rect", xmin = 135, xmax = 197, ymin = 1, ymax = 71.5, fill="white", color="grey50") +  annotate("rect", xmin = 139, xmax = 145, ymin = 9, ymax = 19, alpha = .15, fill="grey50", color="grey50") +  annotate("text", label="Recommended sentence", x=170, y=15, size=3.25) + annotate("text", label="Han defendants", x=162, y=36, size=3.25) +  annotate("text", label="Minority defendants", x=165.5, y=56, size=3.25)+ annotate("text", label="Death sentence", x=12, y=375, size=2.75) + annotate("text", label="Life imprisonment", x=14, y=280, size=2.75) + annotate("text", label="15 Year fixed term imprisonment", x=25, y=195, size=2.75) + geom_point(aes(x=142.5, y=56), colour="brown3") + geom_point(aes(x=142.5, y=36), colour="dodgerblue3")
p2<-ggplot(data.yunnan.heroin.noposs, aes(x=drug.heroin.quantity, y=pun.severity, color=def.minority.fac)) + xlab("Drug Quantity - Heroin (grams)") + ylab("Punishment Severity (months)")  + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE) + theme_bw() + theme(legend.title=element_blank()) + scale_colour_manual(values=c("dodgerblue3","brown3")) + annotate("rect", xmin = 0, xmax = 10, ymin = 0, ymax = 84, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 10, xmax = 50, ymin = 84, ymax = 180, alpha = .15, fill="grey50", color="grey50") + annotate("rect", xmin = 50, xmax = 200, ymin = 180, ymax = 400, alpha = .15, fill="grey50",color="grey50") + geom_point(size=.15) + stat_smooth(method="loess", se=FALSE) + scale_x_continuous(limits = c(0, 200))+ scale_y_continuous(limits = c(-1, 400)) + geom_vline(xintercept = 10, linetype="dashed", size=.8) + geom_vline(xintercept = 50, linetype="dashed", size=.8, color="black") + theme(legend.position = "none") +  annotate("rect", xmin = 135, xmax = 197, ymin = 1, ymax = 71.5, fill="white", color="grey50") +  annotate("rect", xmin = 139, xmax = 145, ymin = 9, ymax = 19, alpha = .15, fill="grey50", color="grey50") +  annotate("text", label="Recommended sentence", x=170, y=15, size=3.25) + annotate("text", label="Han defendants", x=162, y=36, size=3.25) +  annotate("text", label="Minority defendants", x=165.5, y=56, size=3.25)+ annotate("text", label="Death sentence", x=12, y=375, size=2.75) + annotate("text", label="Life imprisonment", x=14, y=280, size=2.75) + annotate("text", label="15 Year fixed term imprisonment", x=25, y=195, size=2.75) + geom_point(aes(x=142.5, y=56), colour="brown3") + geom_point(aes(x=142.5, y=36), colour="dodgerblue3")
arrange_ggplot2(p1,p2, nrow=2, ncol=1)
dev.off() 


##ALL PROVINCE ANALYSIS##

#Table: Estimating Effect of Minority Status 

knots.heroin<-c(10,50)
knots.meth<-c(10,50)
knots.opium<-c(200,1000)
knots.other<-c(10,50)
knots.marijuana<-c(10000,50000) 

m1<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.filter)
m2<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.filter) 
m3<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.filter) 
m4<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.filter) 
m5<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.filter) 
m6<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.filter) 
m7<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.filter) 
m8<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.filter) 

###########################################################################################
#Table A3: Effects of Minority Status on Sentencing: All provinces. "Full Country" columns
###########################################################################################

estimates.pooled.pun.severity <- matrix(data = NA, nrow=8, ncol=7)
colnames(estimates.pooled.pun.severity)<-c("modelnum","estimate", "se", "tstat","pvalue","outcome","covariates")
estimates.pooled.pun.severity<-data.frame(estimates.pooled.pun.severity)
estimates.pooled.pun.severity$modelnum<-seq(1:8)
estimates.pooled.pun.severity$outcome<-rep(c("pun.severity"), c(8))
estimates.pooled.pun.severity$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))

for (j in 1:8) {
  try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
  try(eval(parse(text=paste("t",j,"<-subset(t",j,",t",j,"$term=='def.minority')",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.severity$estimate[",j,"]<-t",j,"[1,2]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.severity$std.error[",j,"]<-t",j,"[1,3]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.severity$pvalue[",j,"]<-t",j,"[1,5]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.severity$n[",j,"]<-df.residual(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(t",j,")",sep=""))))
} 

estimates.pooled.pun.severity$estimate<-as.numeric(estimates.pooled.pun.severity$estimate)
estimates.pooled.pun.severity$se<-as.numeric(estimates.pooled.pun.severity$se)
estimates.pooled.pun.severity$l95ci<-estimates.pooled.pun.severity$estimate-1.96*estimates.pooled.pun.severity$se
estimates.pooled.pun.severity$u95ci<-estimates.pooled.pun.severity$estimate+1.96*estimates.pooled.pun.severity$se


estimates.pooled.pun.severity <- apply(estimates.pooled.pun.severity,2,as.character)
write.csv(estimates.pooled.pun.severity, "estimates.pooled.pun.severity.csv")  

m1<-felm(pun.lifedeath ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.filter)
m2<-felm(pun.lifedeath ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.filter) 
m3<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.filter) 
m4<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.filter) 
m5<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.filter) 
m6<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.filter) 
m7<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.filter) 
m8<-felm(pun.lifedeath ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.filter) 

estimates.pooled.pun.lifedeath <- matrix(data = NA, nrow=8, ncol=7)
colnames(estimates.pooled.pun.lifedeath)<-c("modelnum","estimate", "se", "tstat","pvalue","outcome","covariates")
estimates.pooled.pun.lifedeath<-data.frame(estimates.pooled.pun.lifedeath)
estimates.pooled.pun.lifedeath$modelnum<-seq(1:8)
estimates.pooled.pun.lifedeath$outcome<-rep(c("pun.lifedeath"), c(8))
estimates.pooled.pun.lifedeath$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))

for (j in 1:8) {
  try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
  try(eval(parse(text=paste("t",j,"<-subset(t",j,",t",j,"$term=='def.minority')",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.lifedeath$estimate[",j,"]<-t",j,"[1,2]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.lifedeath$std.error[",j,"]<-t",j,"[1,3]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.lifedeath$pvalue[",j,"]<-t",j,"[1,5]",sep=""))))
  try(eval(parse(text=paste("estimates.pooled.pun.lifedeath$n[",j,"]<-df.residual(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(t",j,")",sep=""))))
} 

estimates.pooled.pun.lifedeath$estimate<-as.numeric(estimates.pooled.pun.lifedeath$estimate)
estimates.pooled.pun.lifedeath$se<-as.numeric(estimates.pooled.pun.lifedeath$se)
estimates.pooled.pun.lifedeath$l95ci<-estimates.pooled.pun.lifedeath$estimate-1.96*estimates.pooled.pun.lifedeath$se
estimates.pooled.pun.lifedeath$u95ci<-estimates.pooled.pun.lifedeath$estimate+1.96*estimates.pooled.pun.lifedeath$se

estimates.pooled.pun.lifedeath <- apply(estimates.pooled.pun.lifedeath,2,as.character)
write.csv(estimates.pooled.pun.lifedeath, "estimates.pooled.pun.lifedeath.csv") 

#Figure: Examining Variation Across Provinces

provinces=list("anhui","beijing","chongqing","fujian","gansu","guangdong","guangxi","guizhou","hainan","hebei","heilongjiang","henan","hubei","hunan","innermongolia","jiangsu","jiangxi","jilin","liaoning","ningxia","qinghai","shaanxi","shandong","shanghai","shanxi","sichuan","tianjin","tibet","xinjiang","yunnan","zhejiang")

for (i in provinces) {
  try(eval(parse(text=paste("m1<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m2<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m3<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m4<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m5<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m6<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m7<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("m8<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.",i,")",sep=""))))
  try(eval(parse(text=paste("estimates.",i," <- matrix(data = NA, nrow=8, ncol=7)",sep=""))))
  try(eval(parse(text=paste("colnames(estimates.",i,")<-c('modelnum','estimate', 'se', 'tstat','pvalue','drug','covariates')",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"<-data.frame(estimates.",i,")",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$modelnum<-seq(1:8)",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$drug<-rep(c('Pooled'), c(8))",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$province<-rep(c('",i,"'), c(8))",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))",sep=""))))
  for (j in 1:8) {
    try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
    try(eval(parse(text=paste("t",j,"<-subset(t",j,",t",j,"$term=='def.minority')",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$estimate[",j,"]<-t",j,"[1,2]",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$se[",j,"]<-t",j,"[1,3]",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$pvalue[",j,"]<-t",j,"[1,5]",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$n[",j,"]<-df.residual(m",j,")",sep=""))))
    try(eval(parse(text=paste("rm(m",j,")",sep=""))))
    try(eval(parse(text=paste("rm(t",j,")",sep=""))))
  } 
  try(eval(parse(text=paste("estimates.",i,"<-data.frame(estimates.",i,")",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$estimate<-as.numeric(paste(estimates.",i,"$estimate))",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$se<-as.numeric(paste(estimates.",i,"$se))",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$pvalue<-as.numeric(paste(estimates.",i,"$pvalue))",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$l95ci<-estimates.",i,"$estimate-1.96*estimates.",i,"$se",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$u95ci<-estimates.",i,"$estimate+1.96*estimates.",i,"$se",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$pvalue.onesided<-estimates.",i,"$pvalue/2",sep=""))))
  try(eval(parse(text=paste("estimates.",i,"$pvalue.onesided[(estimates.",i,"$estimate<0)==TRUE]<-(1-estimates.",i,"$pvalue/2)",sep=""))))
  try(eval(parse(text=paste("estimates.",i," <- apply(estimates.",i,",2,as.character)",sep=""))))
  try(eval(parse(text=paste("write.csv(estimates.",i,", 'estimates.",i,".csv')",sep=""))))
} 

estimates.provinces<-data.frame(rbind(estimates.anhui,estimates.beijing,estimates.chongqing,estimates.fujian,estimates.gansu,estimates.guangdong,estimates.guangxi,estimates.guizhou,estimates.hainan,estimates.hebei,estimates.heilongjiang,estimates.henan,estimates.hubei,estimates.hunan,estimates.innermongolia,estimates.jiangsu,estimates.jiangxi,estimates.jilin,estimates.liaoning,estimates.ningxia,estimates.qinghai,estimates.shaanxi,estimates.shandong,estimates.shanghai,estimates.shanxi,estimates.sichuan,estimates.tianjin,estimates.tibet,estimates.xinjiang,estimates.yunnan,estimates.zhejiang))
estimates.provinces$province.label<-rep(c("Anhui","Beijing","Chongqing","Fujian","Gansu","Guangdong","Guangxi","Guizhou","Hainan","Hebei","Heilongjiang","Henan","Hubei","Hunan","Inner Mongolia","Jiangsu","Jiangxi","Jilin","Liaoning","Ningxia","Qinghai","Shaanxi","Shandong","Shanghai","Shanxi","Sichuan","Tianjin","Tibet","Xinjiang","Yunnan","Zhejiang"), c(9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9))
estimates.provinces$province.label<-ordered(estimates.provinces$province.label, levels = c("Anhui","Beijing","Chongqing","Fujian","Gansu","Guangdong","Guangxi","Guizhou","Hainan","Hebei","Heilongjiang","Henan","Hubei","Hunan","Inner Mongolia","Jiangsu","Jiangxi","Jilin","Liaoning","Ningxia","Qinghai","Shaanxi","Shandong","Shanghai","Shanxi","Sichuan","Tianjin","Tibet","Xinjiang","Yunnan","Zhejiang"))

estimates.provinces$estimate<-as.numeric(paste(estimates.provinces$estimate))
estimates.provinces$se<-as.numeric(paste(estimates.provinces$se))
estimates.provinces$pvalue<-as.numeric(paste(estimates.provinces$pvalue))
estimates.provinces$pvalue.onesided<-as.numeric(paste(estimates.provinces$pvalue.onesided))
estimates.provinces$pvalue.onesided.sig<-"No"
estimates.provinces$pvalue.onesided.sig[estimates.provinces$pvalue.onesided<.05]<-"Yes"
estimates.provinces$l95ci<-estimates.provinces$estimate-1.96*estimates.provinces$se
estimates.provinces$u95ci<-estimates.provinces$estimate+1.96*estimates.provinces$se

write.csv(estimates.provinces, "estimates.provinces.csv")

estimates.provinces<-subset(estimates.provinces, estimates.provinces$modelnum!=7)
estimates.provinces<-subset(estimates.provinces, estimates.provinces$modelnum!=8)

for (i in provinces) {
  try(eval(parse(text=paste("estimates.",i,"<-subset(estimates.provinces, estimates.provinces$province=='",i,"')",sep=""))))
}  

cols <- c("Yes" = "dodgerblue3", "No" = "grey60")
plot.anhui<-ggplot(estimates.anhui, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Anhui", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50') + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=10), legend.text=element_text(size=8)) + theme(legend.position='none')  + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.beijing<-ggplot(estimates.beijing, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Beijing", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.chongqing<-ggplot(estimates.chongqing, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Chongqing", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.fujian<-ggplot(estimates.fujian, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Fujian", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.gansu<-ggplot(estimates.gansu, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Gansu", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.guangdong<-ggplot(estimates.guangdong, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Guangdong", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)  
plot.guangxi<-ggplot(estimates.guangxi, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Guangxi", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + guides(color = FALSE) + theme(legend.position='bottom', legend.box = "horizontal", legend.text=element_text(size=9.75)) + theme(legend.title = element_blank()) 
plot.guizhou<-ggplot(estimates.guizhou, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Guizhou", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))  
plot.hainan<-ggplot(estimates.hainan, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hainan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.hebei<-ggplot(estimates.hebei, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hebei", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.heilongjiang<-ggplot(estimates.heilongjiang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Heilongjiang", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.henan<-ggplot(estimates.henan, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Henan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.hubei<-ggplot(estimates.hubei, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hubei", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.hunan<-ggplot(estimates.hunan, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hunan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.innermongolia<-ggplot(estimates.innermongolia, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Inner Mongolia", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.jiangsu<-ggplot(estimates.jiangsu, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Jiangsu", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.jiangxi<-ggplot(estimates.jiangxi, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Jiangxi", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.jilin<-ggplot(estimates.jilin, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Jilin", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.liaoning<-ggplot(estimates.liaoning, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Liaoning", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.ningxia<-ggplot(estimates.ningxia, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Ningxia", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.qinghai<-ggplot(estimates.qinghai, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Qinghai", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.shaanxi<-ggplot(estimates.shaanxi, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Shaanxi", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.shandong<-ggplot(estimates.shandong, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Shandong", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.shanghai<-ggplot(estimates.shanghai, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Shanghai", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.shanxi<-ggplot(estimates.shanxi, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Shanxi", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.sichuan<-ggplot(estimates.sichuan, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Sichuan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)  
plot.tianjin<-ggplot(estimates.tianjin, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Tianjin", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.tibet<-ggplot(estimates.tibet, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Tibet", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.xinjiang<-ggplot(estimates.xinjiang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Xinjiang", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.yunnan<-ggplot(estimates.yunnan, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Yunnan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2)) 
plot.zhejiang<-ggplot(estimates.zhejiang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Zhejiang", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 

mylegend<-g_legend(plot.guangxi)

#############################################
#Fig A6-Exploring Variations across provinces
#############################################
pdf('fig-provinces-robustness1.pdf', width=9.5, height=10.5)
plot <- grid.arrange(arrangeGrob(plot.anhui + theme(legend.position="none"),plot.beijing + theme(legend.position="none"), plot.chongqing + theme(legend.position="none"),plot.fujian + theme(legend.position="none"), plot.gansu + theme(legend.position="none"),plot.guangdong + theme(legend.position="none"), plot.guangxi + theme(legend.position="none"), plot.guizhou + theme(legend.position="none"), plot.hainan + theme(legend.position="none"),plot.hebei + theme(legend.position="none"), plot.heilongjiang + theme(legend.position="none"),plot.henan + theme(legend.position="none"), plot.hubei + theme(legend.position="none"),plot.hunan + theme(legend.position="none"), plot.innermongolia + theme(legend.position="none"), plot.jiangsu + theme(legend.position="none"),nrow=4,ncol=4),mylegend, nrow=2,heights=c(10, 1))
dev.off()

pdf('fig-provinces-robustness2.pdf', width=9.5, height=10.5)
plot <- grid.arrange(arrangeGrob(plot.jiangxi + theme(legend.position="none"),plot.jilin + theme(legend.position="none"),plot.liaoning + theme(legend.position="none"),plot.ningxia + theme(legend.position="none"),plot.qinghai + theme(legend.position="none"), plot.shaanxi + theme(legend.position="none"),plot.shandong + theme(legend.position="none"), plot.shanghai + theme(legend.position="none"),plot.shanxi + theme(legend.position="none"), plot.sichuan + theme(legend.position="none"), plot.tianjin + theme(legend.position="none"), plot.tibet + theme(legend.position="none"),plot.xinjiang + theme(legend.position="none"),plot.yunnan + theme(legend.position="none"), plot.zhejiang + theme(legend.position="none"),nrow=4,ncol=4),mylegend, nrow=2,heights=c(10, 1))
dev.off()

data.provinces<-read.csv("province.csv",stringsAsFactors=FALSE)
estimates.provinces<-merge(estimates.provinces,data.provinces, by.x="province", by.y="province",all.x =TRUE, all.y=TRUE)
names(estimates.provinces)
names(data.provinces)

####################################################
#Figure A5
####################################################
pdf('fig-estimates-provinces-scatter.pdf', width=6.23, height=6.23) 
ggplot(data.provinces, aes(y=frac.minority, x=log(cases.pc.10000),color=discrimination)) + annotate("rect", xmin = 0, xmax = 2, ymin = 5, ymax = 50, alpha = .15, fill="grey50", color="grey50") + xlim(c(-2.1,2.1))  +  geom_point(size=2.5,alpha = 0) + xlab("Cases per 10000 Citizens (Log)") +ylim(c(0,50)) + ylab("Minority Share of Population (%)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() +  scale_colour_manual(values=c("grey20", "dodgerblue3")) + theme(legend.title=element_blank())+ theme(plot.title = element_text(lineheight=.5)) + geom_text(aes(label=province.label, fontface=bld),hjust=.5, vjust=0, size=3.5) + theme(legend.position = "none")  
dev.off() 

#Figure: Examining Variation Across Groups 

ethnicity<-read.csv("minorities.csv",stringsAsFactors=FALSE)
ethnicity<-ethnicity[1:56,]

data.filter$ind<-1
ethnicity.count<-aggregate(data.filter$ind, by=list(data.filter$def.ethnicity.eng), FUN="sum", na.rm=TRUE)
colnames(ethnicity.count)<-c("minority.english","cases")
ethnicity<-merge(ethnicity,ethnicity.count,by=c("minority.english"), all.x=TRUE,all.y=FALSE) #Consider whether I only want to do filtered cases

ethnicity$pop.2010<-as.numeric(paste(ethnicity$pop.2010))
ethnicity$case.rate<-ethnicity$cases/ethnicity$pop.2010 
drug.meth.quantity<-aggregate(data.filter$drug.meth.quantity, by=list(data.filter$def.ethnicity.eng), FUN="sum", na.rm=TRUE) 
colnames(drug.meth.quantity)<-c("minority.english","drug.meth.quantity")
drug.heroin.quantity<-aggregate(data.filter$drug.heroin.quantity, by=list(data.filter$def.ethnicity.eng), FUN="sum", na.rm=TRUE)
colnames(drug.heroin.quantity)<-c("minority.english","drug.heroin.quantity")

ethnicity<-merge(ethnicity,drug.meth.quantity,by=c("minority.english"), all.x=TRUE,all.y=FALSE)
ethnicity<-merge(ethnicity,drug.heroin.quantity,by=c("minority.english"), all.x=TRUE,all.y=FALSE)

ethnicity$drug.quantity.pooled<-ethnicity$drug.meth.quantity+ethnicity$drug.heroin.quantity
ethnicity$drug.quantity.pc<-ethnicity$drug.quantity.pooled/ethnicity$pop.2010 
ethnicity$drug.share<-ethnicity$drug.quantity.pooled/sum(ethnicity$drug.quantity.pooled,na.rm=TRUE) 

colnames(ethnicity)<-c("minority.english","minority.label","minority.chinese","pop.perc2010","pop.2010","pop.2010.yunnan","pop.2010.yunnan.perc","X","cases","case.rate","drug.meth.quantity.grp","drug.heroin.quantity.grp","drug.quantity.pooled.grp","drug.quantity.pc","drug.share")

data.filter.ethnicity.cc<-subset(data.filter, data.filter$def.ethnicity.eng!="NA" &data.filter$def.ethnicity.eng!="minority")

m1<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.filter.ethnicity.cc)
m2<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.filter.ethnicity.cc) 
m3<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.filter.ethnicity.cc) 
m4<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.filter.ethnicity.cc) 
m5<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.filter.ethnicity.cc) 
m6<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.filter.ethnicity.cc) 
m7<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.filter.ethnicity.cc) 
m8<-felm(pun.severity ~ relevel(as.factor(def.ethnicity.eng),"han") + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.filter.ethnicity.cc) 



for (j in 1:8) {
  try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
} 
t1$label<-gsub(".*)","", t1$term)  
t2$label<-gsub(".*)","", t2$term)  
t3$label<-gsub(".*)","", t3$term)  
t4$label<-gsub(".*)","", t4$term)  
t5$label<-gsub(".*)","", t5$term)  
t6$label<-gsub(".*)","", t6$term)  
t7$label<-gsub(".*)","", t7$term)
t8$label<-gsub(".*)","", t8$term)

minorities=list("achang","bai","blang","bonan","bouyei","choson","dai","daur","deang","derung","dong","dongxiang","ewenki","gaoshan","gelao","gin","han","hani","hezhen","hui","jingpo","jino","kazak","kirgiz","lahu","lhoba","li","lisu","manchu","maonan","miao","monba","mongol","mulao","naxi","nu","oroqen","pumi","qiang","russ","salar","she","sui","tajik","tatar","tu","tujia","uygur","uzbek","wa","xibe","yao","yi","yugur","zang","zhuang")

for (i in minorities) {
  try(eval(parse(text=paste("estimates.",i," <- matrix(data = NA, nrow=8, ncol=8)",sep="")))) 
  try(eval(parse(text=paste("colnames(estimates.",i,")<-c('modelnum','estimate', 'se', 'tstat','pvalue','drug','covariates','outliers')",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"<-data.frame(estimates.",i,")",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"$modelnum<-seq(1:8)",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"$drug<-rep(c('Pooled'), c(8))",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"$province<-rep(c('Pooled'), c(8))",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"$minority.english<-rep(c('",i,"'), c(8))",sep="")))) 
  try(eval(parse(text=paste("estimates.",i,"$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))",sep="")))) 
}
for (i in minorities) {  
  for (j in 1:8) {
    try(eval(parse(text=paste("t",j,".",i,"<-subset(t",j,",t",j,"$label=='",i,"')",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$estimate[",j,"]<-t",j,".",i,"$estimate",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$se[",j,"]<-t",j,".",i,"$std.error",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$pvalue[",j,"]<-t",j,".",i,"$p.value",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$pvalue.onesided<-estimates.",i,"$pvalue/2",sep=""))))
    try(eval(parse(text=paste("estimates.",i,"$pvalue.onesided[(estimates.",i,"$estimate<0)==TRUE]<-(1-estimates.",i,"$pvalue/2)",sep=""))))
  }
}

estimates.minorities<-data.frame(rbind(estimates.achang,estimates.bai,estimates.blang,estimates.bonan,estimates.bouyei,estimates.choson,estimates.dai,estimates.daur,estimates.deang,estimates.derung,estimates.dong,estimates.dongxiang,estimates.ewenki,estimates.gaoshan,estimates.gelao,estimates.gin,estimates.han,estimates.hani,estimates.hezhen,estimates.hui,estimates.jingpo,estimates.jino,estimates.kazak,estimates.kirgiz,estimates.lahu,estimates.lhoba,estimates.li,estimates.lisu,estimates.manchu,estimates.maonan,estimates.miao,estimates.monba,estimates.mongol,estimates.mulao,estimates.naxi,estimates.nu,estimates.oroqen,estimates.pumi,estimates.qiang,estimates.russ,estimates.salar,estimates.she,estimates.sui,estimates.tajik,estimates.tatar,estimates.tu,estimates.tujia,estimates.uygur,estimates.uzbek,estimates.wa,estimates.xibe,estimates.yao,estimates.yi,estimates.yugur,estimates.zang,estimates.zhuang))

estimates.minorities<-merge(estimates.minorities,ethnicity,by=c("minority.english"), all.x=TRUE,all.y=FALSE)
estimates.minorities<-subset(estimates.minorities, is.na(estimates.minorities$estimate)==FALSE & estimates.minorities$pop.2010>400000 & estimates.minorities$cases>25) #Consider increasing this threshold
unique(estimates.minorities$minority.english)

estimates.minorities$estimate<-as.numeric(paste(estimates.minorities$estimate))
estimates.minorities$se<-as.numeric(paste(estimates.minorities$se))
estimates.minorities$pvalue<-as.numeric(paste(estimates.minorities$pvalue))
estimates.minorities$pvalue.onesided<-as.numeric(paste(estimates.minorities$pvalue.onesided))
estimates.minorities$pvalue.onesided.sig<-"No"
estimates.minorities$pvalue.onesided.sig[estimates.minorities$pvalue.onesided<.05]<-"Yes"
estimates.minorities$cases.pc.10000<-estimates.minorities$case.rate*10000

estimates.minorities$l95ci<-estimates.minorities$estimate-1.96*estimates.minorities$se
estimates.minorities$u95ci<-estimates.minorities$estimate+1.96*estimates.minorities$se

write.csv(estimates.minorities, "estimates.minorities.csv")
estimates.minorities<-read.csv("estimates.minorities.csv",stringsAsFactors=FALSE)
estimates.minorities.m5<-subset(estimates.minorities, estimates.minorities$modelnum=="5")

pdf('fig-estimates-minorities-m5.pdf', width=7.23, height=6.25) 
ggplot(estimates.minorities.m5, aes(y=estimate, x=reorder(minority.label,estimate))) + geom_point(size=2.5,alpha = 0.75, color="dodgerblue3") +  geom_errorbar(aes(ymin=l95ci, ymax=u95ci), position=position_dodge(width=.75),width=0, alpha = 0.6, color="dodgerblue3") + xlab("Minority Group") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + coord_flip()+ geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank()) + ggtitle("Outcome: Punishment Severity") + theme(plot.title = element_text(lineheight=.5)) + ylim(c(-22,22))
dev.off() 

pdf('fig-problemminority-full.pdf', width=8.75, height=4.25)  
p1<-ggplot(estimates.minorities.m5, aes(y=estimate, x=log(drug.quantity.pc))) + geom_point(size=2.5,alpha = 0, color="dodgerblue3") + geom_smooth(method=lm,se=FALSE,col="dodgerblue3",alpha=.5) + xlab("Drug Quantity Confiscated per Capita (Log)") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank())+ theme(plot.title = element_text(lineheight=.5)) + geom_text(aes(label=minority.label),hjust=.5, vjust=0, size=3.5)
p2<-ggplot(estimates.minorities.m5, aes(y=estimate, x=log(cases.pc.10000))) + geom_point(size=2.5,alpha = 0, color="dodgerblue3") + geom_smooth(method=lm,se=FALSE,col="dodgerblue3",alpha=.5) + xlab("Cases per 10000 Citizens (Log)") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank())+ theme(plot.title = element_text(lineheight=.5)) + geom_text(aes(label=minority.label),hjust=.5, vjust=0, size=3.5)
arrange_ggplot2(p1,p2, nrow=1, ncol=2)
dev.off() 

pdf('fig-problemminority-full-cases.pdf', width=6.23, height=6.23)  
ggplot(estimates.minorities.m5, aes(y=estimate, x=log(cases.pc.10000))) + geom_point(size=2.5,alpha = 0, color="dodgerblue3") + geom_smooth(method=lm,se=FALSE,col="dodgerblue3",alpha=.5) + xlab("Cases per 10000 Citizens (Log)") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank())+ theme(plot.title = element_text(lineheight=.5)) + geom_text(aes(label=minority.label),hjust=.5, vjust=0, size=3.5)
dev.off() 

#################################################
#Fig2: testing the ``problem minority'' framework
#################################################

pdf('fig-estimates-minorities-final.pdf', width=8.75, height=4.25) 
p1<-ggplot(estimates.minorities.m5, aes(y=estimate, x=reorder(minority.label,estimate))) + geom_point(size=2.5,alpha = 0.75, color="dodgerblue3") +  geom_errorbar(aes(ymin=l95ci, ymax=u95ci), position=position_dodge(width=.75),width=0, alpha = 0.6, color="dodgerblue3") + xlab("Minority Group") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + coord_flip()+ geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank()) + theme(plot.title = element_text(lineheight=.5)) + ylim(c(-22,22))
p2<-ggplot(estimates.minorities.m5, aes(y=estimate, x=log(cases.pc.10000))) + geom_point(size=2.5,alpha = 0, color="dodgerblue3") + geom_smooth(method=lm,se=FALSE,col="dodgerblue3",alpha=.5) + xlab("Cases per 10000 Citizens (Log)") + ylab("Coefficient Estimate on Ethnic Group (months)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + geom_hline(yintercept = 0.0, size=.5, color="grey50")+  scale_colour_manual(values=c("dodgerblue3")) + theme(legend.title=element_blank())+ theme(plot.title = element_text(lineheight=.5)) + geom_text(aes(label=minority.label),hjust=.5, vjust=0, size=3.5) + xlim(c(-2.5,2.5))
arrange_ggplot2(p1,p2, nrow=1, ncol=2)
dev.off() 

summary(lm(estimate~log(drug.quantity.pc),data=estimates.minorities.m5))
summary(lm(estimate~log(cases.pc.10000),data=estimates.minorities.m5))

estimates.minorities<-subset(estimates.minorities, estimates.minorities$modelnum!=7)
estimates.minorities<-subset(estimates.minorities, estimates.minorities$modelnum!=8)

for (i in minorities) {
  try(eval(parse(text=paste("estimates.",i,"<-subset(estimates.minorities, estimates.minorities$minority.english=='",i,"')",sep=""))))
} 

cols <- c("Yes" = "dodgerblue3", "No" = "grey60")
plot.bai<-ggplot(estimates.bai, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Bai", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50') + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=10), legend.text=element_text(size=8)) + theme(legend.position='none')  + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))
plot.bouyei<-ggplot(estimates.bouyei, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Bouyei", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.choson<-ggplot(estimates.choson, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Korean", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.dai<-ggplot(estimates.dai, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Dai", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))
plot.dong<-ggplot(estimates.dong, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Dong", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.dongxiang<-ggplot(estimates.dongxiang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Dongxiang", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))
plot.gelao<-ggplot(estimates.gelao, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Gelao", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + guides(color = FALSE) + theme(legend.position='bottom', legend.box = "horizontal", legend.text=element_text(size=9.75)) + theme(legend.title = element_blank())
plot.hani<-ggplot(estimates.hani, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hani", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.hui<-ggplot(estimates.hui, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Hui", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))
plot.lahu<-ggplot(estimates.lahu, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Lahu", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.li<-ggplot(estimates.li, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Li", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.lisu<-ggplot(estimates.lisu, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Lisu", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.miao<-ggplot(estimates.miao, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Miao", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.mongol<-ggplot(estimates.mongol, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Mongol", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.she<-ggplot(estimates.she, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("She", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.sui<-ggplot(estimates.sui, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Sui", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.tujia<-ggplot(estimates.tujia, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Tujia", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.uygur<-ggplot(estimates.uygur, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Uyghur", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)
plot.wa<-ggplot(estimates.wa, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Wa", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))  
plot.yao<-ggplot(estimates.yao, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Yao", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.yi<-ggplot(estimates.yi, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Yi", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) + theme(panel.border = element_rect(colour = "dodgerblue3", size=2))
plot.zang<-ggplot(estimates.zang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Tibetan", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols) 
plot.zhuang<-ggplot(estimates.zhuang, aes(x=pvalue.onesided, y=estimate, shape=covariates, color=pvalue.onesided.sig)) + xlab('p-value') + ylab('Coefficient Estimate') + theme(plot.title = element_text(lineheight=.6, face='bold')) + theme_bw() + geom_point(alpha=.85, size=2.5) + ggtitle("Zhuang", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color='grey50')  + geom_vline(xintercept =.01,linetype =3,color='grey50') + xlim(0, 1) + ylim(-30, 30) + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position='none') + geom_hline(yintercept =.0,linetype =2,color='grey50') + scale_color_manual(values = cols)

mylegend<-g_legend(plot.guangxi)

#####################################################
#Fig A3 and A4 Robustness of Minority Group Estimates
#####################################################

pdf('fig-minorities-robustness1.pdf', width=9.5, height=10.5)
plot <- grid.arrange(arrangeGrob(plot.bai + theme(legend.position="none"),plot.bouyei + theme(legend.position="none"), plot.choson + theme(legend.position="none"),plot.dai + theme(legend.position="none"), plot.dong + theme(legend.position="none"),plot.dongxiang + theme(legend.position="none"), plot.gelao + theme(legend.position="none"), plot.hani + theme(legend.position="none"), plot.hui + theme(legend.position="none"),plot.lahu + theme(legend.position="none"), plot.li + theme(legend.position="none"),plot.lisu + theme(legend.position="none"), plot.miao + theme(legend.position="none"),plot.mongol + theme(legend.position="none"), plot.she + theme(legend.position="none"), plot.sui + theme(legend.position="none"),nrow=4,ncol=4),mylegend, nrow=2,heights=c(10, 1))
dev.off()

pdf('fig-minorities-robustness2.pdf', width=9.5, height=10.5)
plot <- grid.arrange(arrangeGrob(plot.tujia + theme(legend.position="none"),plot.uygur + theme(legend.position="none"),plot.wa + theme(legend.position="none"),plot.yao + theme(legend.position="none"),plot.yi + theme(legend.position="none"), plot.zang + theme(legend.position="none"),plot.zhuang + theme(legend.position="none"),nrow=4,ncol=4),mylegend, nrow=2,heights=c(10, 1))
dev.off() 

#####################################################
#Table A4: ethnicity.csv
#####################################################
write.csv(ethnicity, "ethnicity.csv")

#NOTE: In Table A4, column (2) "Total Population" corresponds to col "pop.2010" in the csv, col(3) is col(2) divided by 1332810869 (China's total population in 2010, link at http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm), note however, that the total population is a different number according to official government document at: http://www.stats.gov.cn/tjsj/zxfb/201104/t20110428_12705.html/. Last accessed Dec27, 2021.
#col(4) drugs total corresponds to col "drug.quantity.pooled"; col(5) drugs % corresponds to col "drug.share".


##YUNNAN ANALYSIS## 

rm(list=setdiff(ls(), c("data.yunnan","data.yunnan.meth","data.yunnan.heroin","data.court.yunnan.county","arrange_ggplot2","vp.layout")))

data.yunnan$drug.pooled.quantity.sd<-(data.yunnan$drug.pooled.quantity-mean(data.yunnan$drug.pooled.quantity))/sd(data.yunnan$drug.pooled.quantity) 
data.yunnan$drug.pooled.quantity.outlier<-0
data.yunnan$drug.pooled.quantity.outlier[data.yunnan$drug.pooled.quantity.sd>4]<-1
ftable(data.yunnan$drug.pooled.quantity.outlier)

data.yunnan.noout<-subset(data.yunnan, data.yunnan$drug.pooled.quantity.outlier!=1)

#########################################################
#Table 1: Estimating Effect of Minority Status in Yunnan
########################################################
knots.heroin<-c(10,50)
knots.meth<-c(10,50)
knots.opium<-c(200,1000)
knots.other<-c(10,50)
knots.marijuana<-c(10000,50000)

summary(data.yunnan)


m1<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.yunnan)
m2<-felm(pun.severity ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.yunnan) 
m3<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.yunnan) 
m4<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.yunnan) 
m5<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.yunnan) 
m6<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.yunnan) 
m7<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.yunnan) 
m8<-felm(pun.severity ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.yunnan) 


estimates.yunnan.pun.severity <- matrix(data = NA, nrow=8, ncol=7)
colnames(estimates.yunnan.pun.severity)<-c("modelnum","estimate", "se", "tstat","pvalue","outcome","covariates")
estimates.yunnan.pun.severity<-data.frame(estimates.yunnan.pun.severity)
estimates.yunnan.pun.severity$modelnum<-seq(1:8)
estimates.yunnan.pun.severity$outcome<-rep(c("pun.severity"), c(8))
estimates.yunnan.pun.severity$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))

for (j in 1:8) {
  try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
  try(eval(parse(text=paste("t",j,"<-subset(t",j,",t",j,"$term=='def.minority')",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.severity$estimate[",j,"]<-t",j,"[1,2]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.severity$std.error[",j,"]<-t",j,"[1,3]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.severity$pvalue[",j,"]<-t",j,"[1,5]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.severity$n[",j,"]<-df.residual(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(t",j,")",sep=""))))
} 

estimates.yunnan.pun.severity$estimate<-as.numeric(estimates.yunnan.pun.severity$estimate)
estimates.yunnan.pun.severity$se<-as.numeric(estimates.yunnan.pun.severity$std.error)
estimates.yunnan.pun.severity$l95ci<-estimates.yunnan.pun.severity$estimate-1.96*estimates.yunnan.pun.severity$se
estimates.yunnan.pun.severity$u95ci<-estimates.yunnan.pun.severity$estimate+1.96*estimates.yunnan.pun.severity$se

estimates.yunnan.pun.severity <- apply(estimates.yunnan.pun.severity,2,as.character)
write.csv(estimates.yunnan.pun.severity, "estimates.yunnan.pun.severity.csv") 


m1<-felm(pun.lifedeath ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity,data=data.yunnan)
m2<-felm(pun.lifedeath  ~ def.minority + drug.heroin.quantity + drug.meth.quantity + drug.opium.quantity + drug.other.quantity + drug.marijuana.quantity + crime.act + def.recid,data=data.yunnan) 
m3<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid,data=data.yunnan) 
m4<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid | factor(ruling.year),data=data.yunnan) 
m5<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental + def.edu.low | factor(ruling.year),data=data.yunnan) 
m6<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year),data=data.yunnan) 
m7<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(court.name.anon), data=data.yunnan) 
m8<-felm(pun.lifedeath  ~ def.minority + bs(drug.heroin.quantity, knots=knots.heroin,degree=3) + bs(drug.meth.quantity, knots=knots.meth,degree=3) + bs(drug.opium.quantity, knots=knots.opium,degree=3) + bs(drug.marijuana.quantity, knots=knots.marijuana,degree=3) + bs(drug.other.quantity, knots=knots.other,degree=3) + crime.act + def.recid  + def.female + def.age + def.mental +  def.edu.low + def.pleadnotguilty + def.goodattitude | factor(ruling.year) + factor(judge1.name.anon), data=data.yunnan) 


estimates.yunnan.pun.lifedeath <- matrix(data = NA, nrow=8, ncol=7)
colnames(estimates.yunnan.pun.lifedeath)<-c("modelnum","estimate", "se", "tstat","pvalue","outcome","covariates")
estimates.yunnan.pun.lifedeath<-data.frame(estimates.yunnan.pun.lifedeath)
estimates.yunnan.pun.lifedeath$modelnum<-seq(1:8)
estimates.yunnan.pun.lifedeath$outcome<-rep(c("pun.lifedeath"), c(8))
estimates.yunnan.pun.lifedeath$covariates<-rep(c('M1. drug.quantity','M2. M1 + crime.act  + def.recid','M3. M2 + amt cubic spline', 'M4. M3. + Year FE ', 'M5. M4 + def.mental + def.female + def.age + def.edu.low', 'M6. M5 + def.pleadnotguilty + def.goodattitude', 'M7. M6 + Court FE', 'M8. M6 + Judge FE'), c(1,1,1,1,1,1,1,1))

for (j in 1:8) {
  try(eval(parse(text=paste("t",j,"<-tidy(m",j,")",sep=""))))
  try(eval(parse(text=paste("t",j,"<-subset(t",j,",t",j,"$term=='def.minority')",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.lifedeath$estimate[",j,"]<-t",j,"[1,2]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.lifedeath$std.error[",j,"]<-t",j,"[1,3]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.lifedeath$pvalue[",j,"]<-t",j,"[1,5]",sep=""))))
  try(eval(parse(text=paste("estimates.yunnan.pun.lifedeath$n[",j,"]<-df.residual(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(m",j,")",sep=""))))
  try(eval(parse(text=paste("rm(t",j,")",sep=""))))
} 

estimates.yunnan.pun.lifedeath$estimate<-as.numeric(estimates.yunnan.pun.lifedeath$estimate)
estimates.yunnan.pun.lifedeath$se<-as.numeric(estimates.yunnan.pun.lifedeath$se)
estimates.yunnan.pun.lifedeath$l95ci<-estimates.yunnan.pun.lifedeath$estimate-1.96*estimates.yunnan.pun.lifedeath$se
estimates.yunnan.pun.lifedeath$u95ci<-estimates.yunnan.pun.lifedeath$estimate+1.96*estimates.yunnan.pun.lifedeath$se

estimates.yunnan.pun.lifedeath <- apply(estimates.yunnan.pun.lifedeath,2,as.character)
write.csv(estimates.yunnan.pun.lifedeath, "estimates.yunnan.pun.lifedeath.csv") 


